############script for generating latent space DV
rm(list=ls())

library(stringr)
library(igraph)
library(statnet)
library(NetData)
library(dplyr)
library(plyr)
library(tidyverse)
library(hrbrthemes)
library(Cairo)

load("replication/data/before.RData") 
load("replication/data/after.RData")
source("replication/r/gbme.r")
source("replication/r/GetDyadDist.R")


####### start with Pre-ISIS data

### Adjacency matrices used with the GBME function

beforeAdj = data.frame(ego = before$SourceCountry, alter = before$TargetCountry, 
                       valence = before$simon)

beforeAdj = graph.data.frame(beforeAdj, directed=FALSE)
beforeAdj = get.adjacency(beforeAdj, attr="valence")


afterAdj = data.frame(ego = after$SourceCountry, alter = after$TargetCountry, 
                      valence = after$simon)
afterAdj = graph.data.frame(afterAdj, directed=FALSE)
afterAdj = get.adjacency(afterAdj, attr="valence")



##### Load in the GBME output
out=read.table('replication/output/OUT-Be-Sim', header=TRUE)
nss=dim(out)[1]
out=out[round((nrow(out)/10)+1,0):nrow(out), ] # Burn (a burn-in allows for convergence on a value, disregard the burn in period)


# coefficients
coefs <-rbind( apply(out[,-(1:3)],2,mean),apply(out[,-(1:3)],2,sd),
               apply(out[,-(1:3)],2,quantile,probs=c(.025,.975)) )

table = 
  coefs %>% 
  data.frame %>% 
  select(Intercept = b0, 
         Distance = bd1, 
         `Jewish-Pop` = bs1, 
         `Sunni-Pop` = bs2, 
         `Shia-Pop` = bs3, 
         `Polity` = bs4, 
         `GDPpc` = bs5) %>% 
  mutate_all(funs(round(., 2)))

rownames(table) <- c('Posterior Mean', 'SD', '2.5', '97.5')

stargazer::stargazer(table, summary = F, title = 'Coefficient Estimates', 
                     label = 'coef', 
                     out = 'replication/paper/figures/coef.tex')

# convergence diagrams
pDat = 
  out %>% 
  select(Scan = scan, 
         Distance = bd1, 
         `Jewish-Pop` = bs1, 
         `Sunni-Pop` = bs2, 
         `Shia-Pop` = bs3, 
         `Polity` = bs4, 
         `GDPpc` = bs5) %>% 
  gather(key = var, value = val, Distance:`GDPpc`)

# plot diagrams (Figure A.1)
ggplot(pDat, aes(x = Scan, y = val, color = var)) + geom_line() + 
  facet_wrap(~var, scales = 'free_x') + hrbrthemes::theme_ipsum() + 
  labs(x = 'Chains', y = 'Parameter Estimate') + 
  scale_x_continuous(breaks = c(5000, 30000, 50000)) + 
  hrbrthemes::scale_color_ipsum() + 
  theme(legend.position = 'none')
ggsave(filename = 'replication/paper/figures/diagnostics.pdf', 
       device = cairo_pdf)  


# Latent plots
# For a basic latent plot analyis look to 
# example code from Hoff: 
# http://www.stat.washington.edu/hoff/Code/hoff_2005_jasa/gbme.postana.r
z=read.table('replication/output/Z-Be-Sim', header=FALSE)
n<-dim(z)[1]/nss
k<-dim(z)[2]
Pz<-array(dim=c(n,k,nss))
for(i in 1:nss) { Pz[,,i]<-as.matrix(z[ ((i-1)*n+1):(i*n) ,])  }
Pz<-Pz[,,-(1:round(nss/10))]     #drop first 10% for burn in

#find posterior mean of z %*% t(z)
zTz<-matrix(0,n,n)
for(i in 1:dim(Pz)[3] ) { zTz<-zTz+Pz[,,i]%*%t(Pz[,,i]) }
zTz<-zTz/dim(Pz)[3] 
#a configuration that approximates posterior mean of zTz
tmp<-eigen(zTz)
z.pm<-tmp$vec[,1:k]%*%sqrt(diag(tmp$val[1:k])) #use this for distance calculations
#now transform each sample z to a common orientation
for(i in 1:dim(Pz)[3] ) { Pz[,,i]<-proc.rr(Pz[,,i],z.pm) }


rownames(beforeAdj)[18] = 'UAE' # fix UAE name


# plot
z.pm = as.data.frame(z.pm)

pts = data.frame()
for(ii in 1:dim(Pz)[3])
{
  tmp = data.frame(x = Pz[,1,ii],
                        y = Pz[,2,ii], 
                        country = rownames(beforeAdj), 
                        trial = ii)
  pts = rbind(pts, tmp)
}

# geom point
pts$show = ifelse(pts$country %in% c('Saudi Arabia', 'Iran'), 1, 0)
bef_pts = pts



######### Get euclidean distance in before period (this will be used to normalize changes in distance from pre- to post-period producing (FinalData.Rdata))

##get euclidean distance Before
rownames(z.pm) <- rownames(beforeAdj)
BeforeDist <- getDyadDist(z.pm, rownames(z.pm))

#generate dyadic ties
for(j in 1:ncol(BeforeDist)) {
  for(i in 1:ncol(BeforeDist)) {
    before$dyad[i+19*(j-1)] <- BeforeDist[j, i]
  }
}


####get rid of countryA-countryA in edge list 
before <- before[before$SourceCountry != before$TargetCountry, ]

### thin list down to only unique combination of dyads
before <- before[!duplicated(before$dyad), ]

## create scaled dyadic distance variable by diving by average
before$Scaled <- before$dyad/(sum(before$dyad)/nrow(before))


#################


# After ISIS model --------------------------------------------------------



#### load in the GBME output

out=read.table('replication/output/OUT-Af-Sim', header=TRUE)
nss=dim(out)[1]
out=out[round((nrow(out)/10)+1,0):nrow(out), ] # Burn (a burn-in allows for convergence on a value, disregard the burn in period)

# coefficients
coefs <-rbind( apply(out[,-(1:3)],2,mean),apply(out[,-(1:3)],2,sd),
               apply(out[,-(1:3)],2,quantile,probs=c(.025,.975)) )

table = 
  coefs %>% 
  data.frame %>% 
  select(Intercept = b0, 
         Distance = bd1, 
         `Jewish-Pop` = bs1, 
         `Sunni-Pop` = bs2, 
         `Shia-Pop` = bs3, 
         `Polity` = bs4, 
         `GDPpc` = bs5) %>% 
  mutate_all(funs(round(., 2)))

rownames(table) <- c('Posterior Mean', 'SD', '2.5', '97.5')

stargazer::stargazer(table, summary = F, title = 'Coefficient estimates from GBME. Post-ISIL period.', label = 'coef-post', 
                     out = 'replication/paper/figures/coef-post.tex')



z=read.table('replication/output/Z-Af-Sim', header=FALSE)
n<-dim(z)[1]/nss
k<-dim(z)[2]
Pz<-array(dim=c(n,k,nss))
for(i in 1:nss) { Pz[,,i]<-as.matrix(z[ ((i-1)*n+1):(i*n) ,])  }
Pz<-Pz[,,-(1:round(nss/10))]     #drop first 10% for burn in

#find posterior mean of z %*% t(z)
zTz<-matrix(0,n,n)
for(i in 1:dim(Pz)[3] ) { zTz<-zTz+Pz[,,i]%*%t(Pz[,,i]) }
zTz<-zTz/dim(Pz)[3] 
#a configuration that approximates posterior mean of zTz
tmp<-eigen(zTz)
z.pm<-tmp$vec[,1:k]%*%sqrt(diag(tmp$val[1:k])) #use this for distance calculations
#now transform each sample z to a common orientation
for(i in 1:dim(Pz)[3] ) { Pz[,,i]<-proc.rr(Pz[,,i],z.pm) }




##### trying to fix the colors to make top 3 stand out
rownames(beforeAdj)[18] = 'UAE'


# ggplot
pts = data.frame()
for(ii in 1:dim(Pz)[3])
{
  tmp = data.frame(x = Pz[,1,ii],
                   y = Pz[,2,ii], 
                   country = rownames(beforeAdj), 
                   trial = ii)
  pts = rbind(pts, tmp)
}

# geom point
pts$show = ifelse(pts$country %in% c('Saudi Arabia', 'Iran'), 1, 0)

aft_pts = pts

# combine
bef_pts$Period = 'Before ISIL'
aft_pts$Period = 'After ISIL'

pDat = rbind(bef_pts, aft_pts)
pDat$Period = factor(pDat$Period, levels = c('Before ISIL', 'After ISIL'))

# israel and iran
pDat$show = ifelse(pDat$country %in% c('Saudi Arabia', 'Iran'), 1, 0)

# reverse after x-axis for plotting
pDat$x[pDat$Period == 'After ISIL'] = pDat$x[pDat$Period == 'After ISIL']*-1


# FIGURE 2
ggplot(data = pDat, aes(x, y)) + geom_point(color = 'grey', alpha = .2) + 
  hrbrthemes::theme_ipsum() + 
  stat_density2d(data = subset(pDat, show == 1), 
                 aes(x = x, y = y, color = country, 
                     fill = country), 
                 alpha = .2,
                 geom = 'polygon') + 
  scale_color_brewer(type = 'qual', palette = 'Set1') + 
  labs(x = '', y = '') + 
  theme(legend.position = 'bottom', title = element_blank(), 
        legend.text = element_text(size = 12),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.spacing.x = unit(.25, 'cm')) + 
  facet_grid(~Period) + 
  guides(color = guide_legend(nrow = 2, 
                              byrow = T))

ggsave('replication/paper/figures/fig2_iran_saudi-rel-fix.pdf', 
       device = cairo_pdf)


# Figure 3: israel and iran
pDat$show = ifelse(pDat$country %in% c('Israel', 'Iran'), 1, 0)

ggplot(data = pDat, aes(x, y)) + geom_point(color = 'grey', alpha = .2) + 
  hrbrthemes::theme_ipsum() + 
  stat_density2d(data = subset(pDat, show == 1), 
                 aes(x = x, y = y, color = country, 
                     fill = country), 
                 alpha = .2,
                 geom = 'polygon') + 
  scale_color_brewer(type = 'qual', palette = 'Set1') + 
  labs(x = '', y = '') + 
  theme(legend.position = 'bottom', title = element_blank(), 
        legend.text = element_text(size = 12),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.spacing.x = unit(.25, 'cm')) + 
  facet_grid(~Period) + 
  guides(color = guide_legend(nrow = 2, 
                              byrow = T))
ggsave('replication/paper/figures/fig2_iran_israel-REL-fix.pdf', 
       device = cairo_pdf)


pDat$show = ifelse(pDat$country %in% c('Iraq', 'Syria'), 1, 0)


# Figure A.4
ggplot(data = pDat, aes(x, y)) + geom_point(color = 'grey', alpha = .2) + 
  hrbrthemes::theme_ipsum() + 
  stat_density2d(data = subset(pDat, show == 1), 
                 aes(x = x, y = y, color = country, 
                     fill = country), 
                 alpha = .2,
                 geom = 'polygon') + 
  scale_color_brewer(type = 'qual', palette = 'Set1') + 
  labs(x = '', y = '') + 
  theme(legend.position = 'bottom', 
        legend.text = element_text(size = 12),
        title = element_blank(), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.spacing.x = unit(.25, 'cm')) + 
  facet_grid(~Period) + 
  guides(color = guide_legend(nrow = 2, 
                              byrow = T))

ggsave('replication/paper/figures/fig2_iraq_syria-rel-fix.pdf', 
       device = cairo_pdf)


####### Get eulcidean distance in post-IS period

rownames(z.pm) <- rownames(afterAdj)

AfterDist <- getDyadDist(z.pm, rownames(z.pm))

#generate dyadic ties (after)
for(j in 1:ncol(AfterDist)) {
  for(i in 1:ncol(AfterDist)) {
    after$dyad[i+19*(j-1)] <- AfterDist[j, i]
  }
}

### get rid of self-ties
after <- after[after$SourceCountry != after$TargetCountry, ]

### get rid of duplicates in the data
after <- after[!duplicated(after$dyad), ]


### norm dyadic distance variable by dividing by average distance
after$Scaled <- after$dyad/(sum(after$dyad)/nrow(after))

####


#######Generate DV

changeinCoop <- -(after$Scaled - before$Scaled)

#change name and merge
FinalData <- after[ ,c(1,2)]
FinalData$ChangeCoop <- changeinCoop


### save the data as RData
save(FinalData, file = "replication/output/coopDV.RData")
